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Numerically solving the semiconductor Bloch equations within a phenomenological relaxation time 
approximation, we extract both the linear and nonlinear optical conductivities of doped graphene 
and gapped graphene under excitation by a laser pulse. We discuss in detail the dependence of second 
harmonic generation, third harmonic generation, and the Kerr effects on the doping level, the gap, 
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I. INTRODUCTION 

The optical nonlinearity of graphene has been 
predicted^^^ and demonstrated^ to be very strong, which 
makes graphene an exciting new candidate for enhanc¬ 
ing nonlinear optical functionalities in optical devices. 

To optimize the performance of these devices, one of 
the preliminary conditions is to fully understand the 
dependence of the optical nonlinearity of graphene on 
the chemical potential,® temperature, and the excita¬ 
tion frequency. At present, both experiments and the¬ 
ories are still at an early stage. Experiments have in¬ 
vestigated parametric frequency conversion,"^ third har¬ 
monic generation (THG),^®“^® Kerr effects and two- 
photon absorption,second harmonic generation 
(SHG),^®“®^ and two-color coherent control,and 
extracted some third order susceptibilities of graphene 
which are orders of magnitude higher than that of normal 
metal and semiconductor materials. However, the depen¬ 
dence of the nonlinearity on chemical potential, tempera¬ 
ture, and the excitation frequency have not been system¬ 
atically measured. Of the theoretical studies reported, 
most are still at the level of single particle approximations 
within different approaches, which include perturbative 
treatments based on Eermi’s golden rule,®®’^® the quasi- 
classical Boltzmann kinetic approach,and quan¬ 
tum treatments based on semiconductor Bloch equations 
(SBE) or equivalent strategies.When optical tran¬ 
sitions around the Dirac points dominate, analytic ex¬ 
pressions for the third order conductivities can be ob¬ 
tained perturbatively by employing the linear dispersion 
approximation.The calculations show that third 
order conductivities depend strongly on the chemical po¬ 
tential. 

However, there are discrepancies between experimen¬ 
tal results and theoretical predictions. Using the appro¬ 
priate experimental parameters, the susceptibility val¬ 


ues obtained by present theories are orders of magni¬ 
tude smaller than measured values.Possible reasons 
for these discrepancies include: (1) the linear dispersion 
approximation may not be adequate for determining the 
third order nonlinearities; (2) a full band structure cal¬ 
culation beyond the two-band tight-binding model may 
be required; (3) the laser intensity used in experiments 
may be too strong for a perturbative approach, with sat¬ 
uration effects becoming important; (4) thermal effects 
induced by temperature change and gradients may play 
a role in the response, and (5) the inclusion of realistic 
scattering and many-body effects may be required even 
for qualitative agreement with experiment. At a sim¬ 
pler level, different single-particle theories, even based 
on equivalent starting equations at the Dirac cone level, 
have not reached agreement on the final expressions for 
third order conductivities,due to the complexity 
in the analytic calculation. In this work, by numerically 
solving SBE in gapped graphene and doped graphene, we 
address some of these issues by considering the depen¬ 
dence of the optical response on the chemical potential 
and band gap: For weak helds, we investigate whether or 
not the perturbative treatment in our previous work^^ is 
correct and adequate, while for strong fields the numer¬ 
ical results enable us to investigate how saturation can 
affect the nonlinearity. 

We organize this paper as follows: in Sec. H, we present 
our model for a gapped graphene; in Sec. HI, we present 
our numerical scheme in the calculation; in Sec. IV, we 
present our results, which include the comparison to the 
available perturbative formulas and the effects of satura¬ 
tion. We conclude and discuss in Sec. V. 


II. MODEL 

We describe the low energy electronic states by a tight 
binding model, employing pz orbitals 4>a{t',z) with a = 


2 


A, B for different lattice sites. The band Bloch wave 
function of the band can be expanded as 

tpskir, z) = Y^ c“fe$afc(r, z), 


where s is the band index, k = kxX + kyij is the two 
dimensional wave vector, and the Bloch state based on 
site a is 

‘^ak{r-,z) = - Rnni “ Ta,z) . 


Here Rnm = nai + ma 2 is the lattice vector, H is the 
area of one unit cell, ta = 0 and tb = (oi + 0-2) 
are the site positions in one unit cell, and the primitive 

lattice vectors are taken as ai = gq and 

0-2 = cio (^^x + , with the lattice constant qq = 

2.46 A . In our tight binding model, we set the on-site 
energies as A for A sites and — A for B sites, the nearest 
neighbor coupling as 70 = 2.7 eV, and the overlap of the 
Pz orbitals between different sites as zero; the asymmetric 
on-site energies, resulting in a band gap, could be induced 
by a substrate.^® Then the satisfy the Schrbdinger 
equation 


f A 7o/fc\ 

U/fc -aJ 




( 1 ) 


Here /fc = 1 -I- e -I- e is the structure factor. 
The eigen energies and eigenstates are 


Ssk = + (7o|/fc|)^ , s = ±, 

Wfc/ y\/l - A/fe|^ j 



j_ 

•\/2 \ Vl + Afc / 


with A/fc = A/e+k- The band structures for A = 0 and 
0.3 eV are shown in Fig. 1. For nonzero A, the band 
edges are located at the Dirac points K and K', and the 
band gap is 2A. For A = 0, gapped graphene reduces to 
usual graphene, and the low energy dispersion relation is 
massless; for nonzero A the low energy dispersion relation 
is characterized by an effective mass. In the following we 
call A the gap parameter. 

For later use in the discretization of the derivatives in 
Eq. (6), we introduce the matrix elements of as 


J drdzil)l^k^{r,z)e "‘‘"^tljs 2 k 2 {r, z) = 6{ki + q - k 2 ) 

X Us-^ki-,S2ki+q I 


Here Usik-,s 2 k+q is calculated by 

Usik;S 2 k+q = (''siA) '^S2k+q^O‘ik-.a2k+q ■ 



FIG. 1. (Color online) Band structures of gapped graphene 
with different gap parameter A = 0 (red thick curves) and 
0.3 eV (blue thin curves) at A = 0.3 eV. 


with 

Wa.,k-o,2k+q = Y 2 [ drdze-^’t 'T 

^ Rnm '^ol\ : : ■^) ■ 

For small q, we approximate Waik\a 2 k+q ~ 5aia2 

which gives 


Usik-S2k+q — ^ ^ i^sik) ^S2k+q^ 


-tq-Ta 


( 2 ) 


The Berry connections can be found from Us^k\s 2 k+q by 

4siS2fc = * ^qUsik\S2k+q)\q—Q > (3) 

and then the velocity matrix elements are given as Vg^k = 
fr^VkSsk and Usjfc = ih~^{£sk-£sk)^ssk, with s = -l-(-) 
when s = — (+). After some algebra, we find 

v+-k = {c+kT <^-k9k + {c+kT c-k9k 

= ^ {*Ini[/fcgfc] + AfkRelfkgk]} , (4) 

with Qk = fi-^'yo[Vkfk + i{TB - TA)fk]- We are in¬ 
terested in optical transitions around the Dirac points 
K = (6i -I- 262)73 and K' = (26i -I- 62)73 with the prim¬ 
itive reciprocal lattice vectors 61 = ^ ~ y) and 

62 = + The usual approximated quanti¬ 

ties around the Dirac points which we used are listed in 
Table I. 

With the application of an external homogeneous elec¬ 
tric field E{t), within the independent particle approxi¬ 
mation the time evolution of the system can be described 
by SBE37 

-dpk{t) 


ih- 


dt 


= [£k - eE(t) ■ Ifc, pk{t)] - ieE{t) ■ VkPkit) 


+ ih 


dpk{t) 


dt 


( 5 ) 
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TABLE I. Lowest order approximations for k around the 
Dirac points. Here we use Vf = /{2h), k = kk 

with k = cos9x + sin Oy, = hvpi^lyj + {hvFn)^, and 
pK = A/^/A2 + {hVpKy 



k — K + K 

k = K' + K 

'70 fk 

ihvFKe~''^ 

ihvFne’’^ 

£ + k 

'^ + + h 

y/{hVFKY + A 2 , 
anVpk 

yJihVFKY + A 2 

anVpk 

Qk 

VF(ix + y) 

VF{ix — y) 

■Si ■Si 

1 1 

wfc ”®® (i/Sk cos 9 — sin 9) 
wfc ”*® {IPk, sin 9 + cos 9) 

ufc ®® (i^K cos 9 + sin 9) 
VFe*® {iPk, sin 9 — cos 9) 


Here pk is a single particle density matrix, for which the 
diagonal term pssk gives the occupation at state ipsk and 
the off-diagonal term p_|_ k identifies the interband po¬ 

larization between two bands; £k is the energy matrix 
with elements £sis-ik = < 5 siS 2 ^sife; and e = — |e| is the 
electron charge. Although ^k alone is a gauge dependent 
quantity, depending on the phases chosen for the Bloch 
functions, the combination with the derivative term Vfc 
is gauge independent and can be written as 


[-eE{t) ■ Ifc, Pk] - ieE{t) ■ VkPk 

icEi^t') • Vq {Uk-k+qPk+qUk+q;k)\q—Q j (6) 


The term 

at 


describes the relaxation processes. 


In a phenomenological way we model the intraband (in¬ 
terband) relaxation process by a parameter (Fg), and 
then 


^Pssk (0 


dt 

dpssk (0 


dt 


ri[Pssfc(0 Pssfe] ) 
^ePssk {t') , 


(7) 


where the density matrix at equilibrium state is given by 
Psisafc = [7 + at temperature T 

and chemical potential p. The current density is calcu¬ 
lated as 


SlS2 


dk 

(2^)2 


'Vs2SikPsiS2k(t) . 


( 8 ) 


To focus on the nonlinear response, we separate the lin¬ 
ear and nonlinear contributions to the perturbed density 
matrix by writing pk(t) = Pk + Pk^(t) + p^k'^'’(t) where 
Pk'^ (t) is the perturbative linear contribution of the elec¬ 
tric field and determined by 


ih 


dt 


[£k,pi^\t)] - eE{t) ■ {[a,pO] -f iVkpl} 


—i 


VeP^-ikit) 


^eP^^Ut)] 


(9) 


while (t) includes all higher order contributions and 
satisfies the equation 


ih 


dp^:^\t) 

dt 


[£k,p‘C‘\t)] - eE{t) ■ + p‘'k'\t)] 

+i'^ k[Pk\t) + Pk'^Ht)]] 


—i 


(^^Pf+kit) 

Vept^t) 


r^p’L'Mt)) 


( 10 ) 


The solution of Eqs. (9) and (10) completely determines 
the evolution of the single-particle density matrix, and 
the current can be written as J‘^{t) = ji^^’'^(t)-|-J*^"*^’'^(t) 

where and ji”*)’‘^(t) are induced by p^k\^) 

p^k^\t) respectively, and describe the linear and nonlin¬ 
ear response. 


III. NUMERICAL SCHEME AND FITTING 
PROCEDURE 


We consider the response of the current to an applied 
electric field pulse with a Gaussian envelope function, 

E{t) = -F c.c., (11) 


with a duration Ac and a center frequency Wc- In the 
frequency domain, this corresponds to a function with 
Gaussian peaks at iwc 


E{uj) = j dte^^*E{t) 

= y^AcEo 




each with a spectral width 21 Ac. 

In contrast to the numerical study by Zhang et 
where the p ■ A interaction is used and there is no cou¬ 
pling between different k points, our SBE, which is based 
on the r ■ E interaction, involve derivatives of the single¬ 
particle density matrix with respect to fc. In the numeri¬ 
cal calculations, we divide the Brillouin zone (BZ) into an 
M X M homogeneous grid, and discretize the derivative 
in Eq. ( 6 ) as^° 


'^<iFiq)\q=o 


V— 

167r2 ^ M \m) ’ 


( 12 ) 


where Qi are chosen as six symmetric points of the hon¬ 
eycomb lattice 


{^ 1 , b2, —bi, — 62 , bi - 1 - b2, —(hi -I- ^ 2 )} , 


Throughout this work, we are interested in the opti¬ 
cal response at different frequencies and its dependence 
on the electric field amplitude Eq, the chemical poten¬ 
tial p and the gap parameter A. Other parameters used 
in the simulation are fixed as T = 300 K, Ac = 100 fs, 
huc = 0.6 eV, and F^ = Fe = 33 meV. The discrete k 
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points are taken from a grid with M = 1500, and in¬ 
cluded in the calculation if e+k < 3.5 eV; tests involving 
the inclusion of more k points confirm that such a restric¬ 
tion leads to converged numerical simulations. The time 
evolution of Eqs. (9) and (10) is solved by a fourth order 
Runge-Kutta method with a time step At = 0.05 fs. The 
current in Eq. (8) is numerically calculated by summing 
all band indices and all the effective k points on the grid 
with an equal weight. After discretization, Eqs. (9) and 
(10) become linear differential equations for which the 
accuracy of the numerical solution is only limited by the 
time step. We point out that the density matrix pfc(t) ac¬ 
quires a phase dependence on k that changes with time. 
At long enough times pk (t) can be strongly dependent on 
fc, and then an accurate calculation of the current from 
Eq. (8) requires a very dense grid, without which the 
nonlinear current is buried in numerical noise. Similarly, 
a dense grid is also required if the relaxation parameters 
Tj/g are very small. However, when making calculations 
for the pulses and relaxation parameters we adopt here, 
we find that the nonlinear current can be determined re¬ 
liably by the use of the moderate grid identified above. 




huj (eV) 

FIG. 2. (Color online) Linear optical current (a) 

(b) Jbh- (cj). The parameters used in the calculation are 
Eq = 10® V/m, A = 0.10 eV, and p = 0. In (b), squares are 
numerical results, while the curves are fitted to = 

ai{uJc)E'^{Lj). 

We begin by illustrating the fitting procedure used in 
this work to extract the coefficients characterizing the 
optical response, and consider a weak incident optical 
pulse with Eq = 10® V/m, A = 0.10 eV, and ^ = 0. 
The linear response can be determined by solving Eq. (9) 
numerically, and using the result to construct 
The result is shown in Fig. 2(a) for an incident field in 
the X direction, and the Fourier transform. 




hcu (eV) 

FIG. 3. (Color online) Nonlinear optical current (a) 
and (oj) for us around (b) Uc, (c) 2(jJc, (d) 3 cJc. The 

parameters used in the calculation are Eq = 10® V/m, A = 
0.10 eV, and p — 0. In figures (b)-(d), squares are numerical 
results, while the curves are fitted to Eq. (15), and the fitting 
parameters are given in the text. 

is numerically determined and shown in Fig. 2(b). Very 
generally the linear response is of the form 

and could be extracted directly for w around 

ojc- Putting ai{oJc) = the result is ui{iOc.) = 

(1.11 — 0.05f)cro, with the universal conductivity (Tq = 
e^/ (41i). 

The situation is different for the nonlinear response. 
It can be determined by solving Eq. (10) numerically, 
and using the results to construct The result is 

shown in Fig. 3(a); note that it is much smaller than the 
linear response, and the peak amplitude is shifted to a 
time slightly later than the peak of the linear response. 
The Fourier transform, 

can then be numerically determined. Here we find a sig¬ 
nificant response for oj close to Wg [corresponding gener- 
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ally to the Kerr effect and two-photon absorption, and 
shown in Fig. 3(b)], for w close to 2ujc [corresponding to 
SHG, and shown in Fig. 3(c)], and for uj close to 3wc 
[corresponding to THG, and shown in Fig. 3(d)], and 
of course for a; close to the associated negative frequen¬ 
cies. While Figs. 3(c) and 3(d) are essentially Gaussian in 
form, as was Fig. 2(b), Fig. 3(b) certainly is not. So the 
question arises as to how to characterize the nonlinear 
response and identify the relevant response coefficients. 

In the perturbative regime, we can very generally ex¬ 
pect a nonlinear response of the form 

= J - wi) 

xE^{oji)E^{oj2)E^{uj-uJi-oj2). (13) 


For LOij^k = iwc and small an approximate analytic 
perturbation calculation^"^ leads to 


(2):xxx 
(j\ n 

^(3);xxxx 


, „ r \ (2) leV 

(wi -l- (5i, Wj -l- (52 j Wfe -|- da) ~ ^ 


leV 


„(3) 


^■(di -b d2 -b da) -b fy'" 



(14) 


(2) (3) 

where s\ , s) , and 7 are determined in the calculation, 
and take on different values for different choices of the 
uji,j,k- Motivated by this, we fit the results of Fig. 3 by 
assuming the conductivity as the form of Eq. (14) with 
taking s\^\ and 7 as free fitting parameters. This 
leads to a fit of the nonlinear current spectrum around 
fl = toc, 2uJc, and 3wc of the form 


10“^®m^V^, cr(C^S 2 ^^ = (0.9 —202.5z) x 10“^® m^/V^, and 
7 = 36.3 meV. 

For weak incident fields, we can use this strategy to 
extract coefficients Wc), Wc, Wc) 

and from our numerical calcula¬ 

tions, confirm that they are independent of the amplitude 
Eq of the incident field - as they should be in the pertur¬ 
bative regime - and compare them with the results of the 
approximate but analytic expressions for these response 
coefficients. For strong incident fields a strict perturba¬ 
tive response of the form in Eq. (14) is not expected to 
hold. Still, the nonlinear response can be expected to 
be characterized by SHG, THG, and terms that behave 
phenomenologically as Kerr and two-photon absorption 
effects. Thus from our numerical calculations we can 
extract an effective Wc) [which we denote as 

o’shg(<x’c)], an effective Wc, Wc) [which we de¬ 

note as (Tthg(<x’c)] and an effective a^^'>’^^^^{—uJc,uJc,LOc) 
[which we denote as ani{i^c)]- Unlike the coefficients that 
govern the perturbative regime, we can expect the ef¬ 
fective coefficients crsHG(wc), crTHG(wc), and ani{u!c) to 
depend on the amplitude of the electric field strength, 
containing renormalizations of the perturbative response 
coefficients in the presence of strong fields. 

Using the Htting scheme described above, we study two 
examples of the dependence of the effective conductivities 
on the chemical potential /x, the gap parameter A, and 
the electric Held amplitude Eq. In the first we consider 
the dependence on /x and Eq with A = 0, which we refer 
to as doped graphene (DG). In the second we consider 
the dependence on A and Eq with /x = 0, which we refer 
to as undoped gapped graphene (GG). 


IV. RESULTS 


j(");-(L! + d) = Cn 


(n 

Si 


leV 

_ ( 

hS + X 7 




y/n 


(15) 


Here n = 2 is used for H = 2wc, and rx = 3 for U = Wc or 
3wc, with Co describing the permutation factor relevant 
for the nonlinear process; C^^^ = 3, and C2ujc = = 1. 

The result of this fitting is shown by the solid curves 
in Fig. 3(b)-3(d), and we can see that indeed a very 
good fit is provided. Once the fit of Eq. (15) is ac¬ 
cepted, we can return to Eq. (14) and identify the non¬ 
linear response coefficients (associated 

with SHG for a fundamental at Wc), 

(associated with THG for a fundamental at Wc), and 
(T(3);^^'^“(_a;c,a;c,Wc) (associated with the Kerr effect 
and two-photon absorption for a fundamental at ojc)- 
For the results shown in Fig. 3, for example, we find 
aQ^a^^^’^^^(oJc,^^c) = (—37.2 + 16.5i) pm/V, with ~ 
0; aQ^a(^^’^^^^(uJc,uJc,u;c) = (0.91-b0.04f) x IQ-i® mVV^ 
with s® ~ 0; and Wc, Wc, Wc) = (—4.1 -b 

0.2f)x 10-16 mVv2 with = (1463.5-b 239.7x) x 


A. Comparing numerical calculations to analytic 
perturbation results 


As a benchmark, we first compare the numerical effec¬ 
tive conductivities at a weak electric field Eg = 10® V/m 
with those available from analytic perturbation calcula¬ 
tions. We begin with the linear response. For DG, in 
previous work3’^ we presented the analytic expression for 
(t(i);'^“(w) obtained perturbatively from the same SHE as 
Eq. (5), taking into account both interband and intra¬ 
band relaxation coefficients Fg and Fj respectively, but 
using matrix elements and energies correct only around 
the Dirac points; our analytic result is 

/ OO 

Ef,(x, T)ll-E^(x, T)]cr^^l^.g (uj; x)dx , 

-00 

(16) 

where /3 = IKUbT) with ks Boltzmann’s constant, and 
E^(x,T) = [1 -b The conductivity at zero 

temperature is 


'^dgT(^^5 m) = — <{ -Q +\(^ + fBg) -b 


XCTo 


4|mI 


hw + iVi 
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FIG. 4. (Color online) The linear effective conductivities for 
DG (squares) and GG (dots) for Eq = 10® V/m. The curves 
are calculated from Eq. (16) for DG and Eq. (17) for GG. 


Here the function (0) is given for 9 = 9,. + i9i as 


g\f,\i9)= In 


2|^| + 9 


— arctan 


2\^l\-9 
9r T 2|^| 


+ t TT + arctan 


9r-2\^x\ 


For GG, because the chemical potential is taken as 0 
and the gap is nonzero, the net contribution to the linear 
conductivity from the intraband transitions (Drude term) 
vanishes at zero temperature; even at room temperature 
that contribution is negligible, so for GG we can restrict 
the expression for the linear conductivity to its interband 
component. 


a 


{1)\XX 

GG;inter 


(w) 



- ^sk) 

i^Ssk “t“ ^Fg 


icrp 

TT 


+ iV e) + 


4A 


hui + i+e 


icrp 

TT 


(2A)^ 

{huj + iVe)'^ 


+ iFe). 


(17) 


In Fig. 4 we plot the results extracted from our numerical 
simulations of Eq. (9), together with the analytic results 
in Eq. (16) and (17) as a function of /x (for DG) and A 
(for GG). The agreement is very good. 

Turning next to the third order response, for 
the analytic expressions of Wg,Wg,Wg) and 

tT(3);^®^“(a;c, Wg, Wg) relevant for DG we use our previous 
results,including both interband and intraband relax¬ 
ation, and with matrix elements and energies taken to 
be those that characterize the regions about the Dirac 
points. For GG with a nonzero gap parameter, pertur¬ 
bative results for THG were obtained by Jafari,^° but in¬ 
stead of using the SHE in Eq. (5) a Kubo formula based 
on the p ■ A interaction was used, without the inclusion 
of any relaxation. Thus while we present our numerical 


results for cr„/(wc) and tTTHG(^^c) for both DG and GG, 
we only compare with the relevant analytic results from 
perturbation theory obtained for DG. This is shown in 
Figs. 5(a) and 5(b) for Gnii+c) and crTHG(a^c) respectively 
at Ep = 10® V/m. The numerical and analytic results 
for DG match very well for chemical potentials over the 
range shown. There is a noticeable difference between 
the numerical and analytically results for Re[cr„i(a;g)], 
although it is less than 10%, for p < 0.3 eV. We can 
attribute this to the singular behavior that Re[cr„i(ajg)] 
exhibits in the perturbative calculation®^ for |/i| < &Ug/2, 
here \p\ < 0.3 eV. Associated with this, the nonlinear 
current in the numerical calculation shows a very strong 
dependence on the pulse duration and shape, and the 
strategy identified above for extracting ani(uJc) from the 
pulse calculation is not completely successful. 

The very good agreement at Ep = 10® V/m between 
the effective conductivities of DG extracted from the nu¬ 
merical calculations, and the conductivities predicted by 
the analytic perturbation theory, suggests that Eqs. (14) 
and (15) provide a reasonable fitting procedure, and as 
well that for weak fields the perturbative results pre¬ 
sented earlier®^ are reliable. It also indicates that the 
usual Dirac point approximations adopted in the pertur¬ 
bative calculation, involving the linear dispersion relation 
and the form of the matrix elements, do not introduce 
any significant errors in calculating the linear and non¬ 
linear optical response of DG at incident photon energies 
around hw = 0.6 eV. 

We also see from Figs. 5(a) and 5(b) that there is a 
similarity in the dependence of the DG results on p with 
the dependence of the GG results on A. Before turning 
to the response of both system at larger field strengths, 
we address such similarities in the following section. 


B. Comparing DG and GG 

In investigations of the optical conductivities of doped 
graphene, 2\p\ is often treated as an effective gap.®’®^ 
Since GG has a real gap of 2A, it is interesting to com¬ 
pare the dependence of the optical conductivities on the 
effective gap 2\p\ induced by the chemical potential in 
DG with the real gap 2A arising in GG. In linear re¬ 
sponse, some insight can be gleaned by comparing the 
analytic formulas in Eq. (16) and (17) for DG and GG. 

In Eq. (17) the interband velocity matrix elements t>-|_ k 

depend on (3^, as shown in Table I, and through that de¬ 
pendence they depend on A. If (3^. were not present, only 
the hrst term in the bracket of Eq. (17) would survive, 
corresponding to the interband contribution to the con¬ 
ductivity of DG®^ with |/i| replaced by A. The presence 
of (3k. leads to the appearance of the other two contribu¬ 
tions. Interestingly, one has the same form as the Drude 
term in DG (with \p\ replaced by A), while the other is 
new. 

The consequences of the second new term in GG are 
apparent in the results shown in Eig. 4; the main dif- 
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FIG. 5. (Color online) Nonlinear response for DG and GG: the nonlinear conductivity a-Q^aniii^c) and cr^^crTHG(t^c) at (a), 
(b) Eo = 10® V/m and (e), (f) 2 x 10^ V/m; of GG with A = 0.20, 0.25, 0.30, and 0.35 eV for uj around (c) u)c and 

(d) 3ujc. The y-axis for the real (imaginary) parts of S'f the left (right) hand side of (c) and (d). Solid curves are 

calculated from analytic perturbation results®^ for DG; dashed curves are drawn to guide the eye. 


ferences between the DG and GG results is that around 
huJc ~ 2A the latter show a deeper valley in the imag¬ 
inary part of the conductivity, and a larger peak in the 
real part of the conductivity. The real part of the conduc¬ 
tivity is associated with absorption, and through Fermi’s 
Golden Rule it is determined by both the joint density 
of states and the velocity matrix elements. Now for GG 
the joint density of states for energies around the Dirac 
points is given by 

= 2 X! = 2-K(hvFf^^^ ~ ’ 

where the factor 2 comes from the spin degeneracy. Gom- 
paring with DG, the joint density of states is the same for 
e > 2|/r| in DG as it is for e > 2A in GG, and so at such 
energies the differences in the linear conductivity should 
be associated with the velocity matrix elements; and in¬ 
deed, they can be linked to the last term in Eq. (17). 

Turning to the nonlinear response, we first compare the 
DG and GG results for cr„;(wc), shown in Fig. 5(a). The 
result for GG shows hne structure as A is close to tujjcl2. 
For A < huJc/2 both one- and two-photon absorption are 
present, and Re[CT„i(wc)] is negative and increases in mag¬ 
nitude with increasing A. In a manner similar to what 
is shown by the results of perturbative calculations^’^^ at 
\^\ < hwcl2 for DG, we expect that at A < fia;c/2 for GG 
the two-photon absorption is associated with saturation 
as described at the level of the third-order nonlinearity; 


it would diverge when relaxation effects are not included. 
For A > huJc/2^ where only two-photon absorption exists, 
the negative value of Re[CT„i(wc)] is induced by the inclu¬ 
sion of the relaxation.Maximum absolute values of the 
imaginary and real parts of ani(uJc) occur for GG around 
A = fu^cl2, and the differences between the results for 
GG and DG can again be attributed to the velocity ma¬ 
trix elements. 

We turn to the results for crTHG(wc) shown in Fig. 5(b). 
The expected similarity of the results for GG and DG, 
respectively as a function of A and |/i|, fails mainly for 
A,^ > 0.25 eV. Here Re [tTTHG(wc)] for GG increases 
faster than that of DG, as functions of A and |^| respec¬ 
tively, while the dependences of Im[(TTHG(<4’c)] for GG 
and DG are analogous, but with larger absolute values 
for GG. Again these differences can be traced back to the 
different velocity matrix elements. 

Here we shortly discuss the relation between the fit¬ 
ted effective conductivity at Wc and the amplitude of 
the optical current calculated from a laser pulse. As in 
Eq. (14), the conductivity shows a strong frequency de¬ 
pendence, and thus the value of the conductivity ani{ujc) 
at the center frequency of a light pulse is generally not a 
good indication of the amplitude of the optical response 
j(nO;x(^) exciting pulse of light is actually used. 

The numerical results for are shown for GG 

at (jj close to Wc in Fig. 5(c), and for w close to Swc 
in Fig. 5(d), for A = 0.20, 0.25, 0.30, and 0.35 eV. At 



















A = fiwc/S = 0.30 eV, both the real and imaginary parts 
of the nonlinear optical current [black curves in Fig. 5(c)] 
show a different shape than those at other A, although 
they do not really exceed them in amplitude. In con¬ 
trast, there is no obvious shape distortion in the spectrum 
shown in Fig. 5(d). As such, the values of crTHG(wc) are 
consistent with the magnitude of the optical current of 
the THG components. 

The results for ani{u!c) and crTHG(wc) extracted from 
the numerical results for a larger Eq = 2 x 10^ V/m 
are shown for both GG and DG in Fig. 5(e) and 5(f). 
Note that the dependence of the effective coefficients of 
DG on \fi\, and those of GG on A, are similar in na¬ 
ture to the dependence of those effective coefficients at 
Eq — 10® V/m, but they take on different values. Hence 
we are now beyond the perturbative regime, and can¬ 
not link the effective coefficients ani{uic) and crTHG(wc) 
with the perturbative results for Wc,Wc,Wc) 

and respectively. For ani{<^c) there 

are significant differences between the values at Eq = 
10® V/m and at Aq = 2 x 10^ V/m at all values of |/i| 
(or A), while for crTHG(wc) the differences are substantial 
only for A, |/r| < 0.3 eV. We attribute these differences to 
saturation effects, which we discuss in the next section. 


C. Saturation effects 


We now turn to the dependence of the effective co¬ 
efficients ani{uJc) and crTHG(<j-'c) on field strength. We 
begin with (Tni(ujc), and note that there are two different 
regimes that we can identify for both DG and GG: 

(i) 2|^| < HlUc for DG, or 2A < huic for GG. Here 
one photon absorption exists and carriers can be injected 
from the ” band to the “-I-” band. Stronger electric 
fields inject more carriers. If the electrons have a finite 
lifetime in the states into which they are injected, their 
injection prevents the effectiveness of further absorption. 
Phenomenologically, the effect of the injected carriers on 
the total absorption a is often characterized by introduc¬ 
ing a saturation field strength Agat, 


OtQ 


1 + {E/Es,tY 


(18) 


with Go the linear absorption and E the electric field am¬ 
plitude in an assumed continuous wave excitation at fre¬ 
quency oj. For isolated graphene, the absorption of nor¬ 
mally incident light is proportional to Re[cr|g(a;)], where 
cr/f (w) is a field dependent effective conductivity, and we 
would expect 


ReKKo;)] = 


Re[( 


(l)-,XX 




\2 ' 


(19) 


1 + (A/Agat) 

However, at weak fields we have^^ 

<^(u;) =a(^)’^^(w)-b3CT(®)’““““(a;,a;,-a;)F;2, (20) 


where w, —w) is the third order conductiv¬ 

ity resulting from a perturbative calculation. Comparing 
with the weak field expansion of Eq. (19) we find 


Ega t — \ I 


Re[cr(i)’““(a;)] 


3Re[(T(®h^^^®(—w, u), w)] 


( 21 ) 


For strong electric fields, we assume that Eqs. (19) and 
(20) work for a field dependent conductivity ani{u!c)', fur¬ 
ther, since we extract ani{ujc) from a numerical calcula¬ 
tion with the incident field in Eq. (11) we can identify 


for Ac > hjVi^e where the pulsed excitation approaches 
continuous wave excitation, and we can replace E by Eq 
in Eq. (19); then we find 


Re[(T„/(w)] 


Re[^(l);a:x(^)] j 

1 + (Eo/Egat)' ■ 


( 22 ) 
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FIG. 6. (Color online) Electric field dependence of the non¬ 
linear conductivities Oniiujc) for different (p. A), (a) (0,0) 
(squares), (0.2,0) eV (circles), and (0,0.2) eV (diamonds), 
(b) (0.4, 0) (squares) and (0, 0.4) eV (diamonds). The real and 
imaginary parts are given by filled and hollow symbols, and 
their scales are at the left and right y axis, respectively. The 
solid curves in (a) are fitted by functions — ^ 

with two fitting parameters cr/ and Egat, while the dashed 
curves are drawn to guide the eye. 


In Fig. 6(a) we plot the dependence of the numeri¬ 
cally determined (Jni(oJc) as a function of Eq for three 
different parameter sets (fi, A) = (0,0), (0.2,0), and 
(0,0.2) eV. The real part of (t„/(wc) is fitted to an ex¬ 
pression —o-//[3(if/at + Eg)] with two parameters cr/ and 
Egat ■ We find that the fittings (shown as solid curves) are 
very good for Re[cr/]/(To ~ —1, —1, and —1.53 respec¬ 
tively, with the saturation fields the same in all cases as 
Egat ~ 3 X 10^ V/m. Gomparing the fitted form with 
Eq. (22), and noting that the linear conductivities are 










9 


given by cr(^)’““(wc)/o'o = 1, 0.96 —0.Hi, and 1.38 —0.36i 
respectively for our parameter sets, the closeness of the 
fitted tJ/ with these linear conductivity indi¬ 

cates that the saturation can indeed be attributed to lin¬ 
ear absorption. Further, by using the numerical values of 
Cni (wc )/cto at a weak field value Eq = 10® V/m, which for 
our three parameter sets are —3.7 x 10“^®, —3.7 x 10“^®, 
and —5.7 x 10“^® m^/V^ respectively, Eq. (21) leads to 
saturation fields of 3 x 10^, 3 x 10^, and 2.8 x 10^ V/m, 
which are very close to the fitted values. The field de¬ 
pendence of Im[(T„/(wc)], which at least in the weak field 
limit can be related to the real part of the nonlinear re¬ 
sponse via nonlinear Kramers-Kronig relations, varies in 
a more complicated way. 

The saturation field can also be estimated from only 
the linear absorption coefficients. Physically, the satu¬ 
ration effect occurs when the injected electron density 
from one-photon absorption is comparable to the density 
of states in the region of k space where the electrons are 
injected. The injected electron density is h/Ti^^^{ujc)E'^ 
with the one-photon absorption coefficients® = 

2 Re[(T(^^’'^^(a;c)]/(?ia;c) and the critical field amplitude 
Em, while the total available states are estimated as those 
satisfying —Tg < s+k — £-k — kuc < Tg, which has a den- 
Then the critical field amplitude 

Em is estimated as 


E„ 


l2T,r, 


o-Q 


huJr 


TT Re[cr(^k^^(a;c)] h\e\vF 


(23) 


This can be used to hnd approximate values of Em 
2.8 X 10^ V/m for those two parameter sets considered 
for DG, and Em ^ 2.4 x 10^ V/m for the parameter set 
considered for GG. Both values are close to the fitted 
saturation field. 

(ii) 2\fi\ > hujc for DG, or 2A > hwc for GG. Here 
we focus on the frequency regimes 2|^| > fiwc > |/r| 
or 2 A > HllIc > A where two photon absorption ex¬ 
ists. Two photon absorption can inject carriers, but it 
is less efficient than one photon absorption. Thus satu¬ 
ration requires higher electric fields, and Eq. (19) does 
not correctly describe the physics, as shown in Fig. 6(b) 
for two parameter sets = (0.4,0) and (0,0.4) eV, 

which has different tendencies compared to the curves in 
Fig. 6(a). For the electric field up to Eq = 2 x 10® V/m, 
the imaginary part of tT„;(wc) does not change much for 
either of these examples. The real part of CT„i(wc) of DG 
changes from negative values to positive values around 
Eq ^ 4: X 10 ^ V/m; while that of GG remains positive 
and decreases. For photon energies where even two pho¬ 
ton absorption is absent, we believe that saturation can 
only occur for much higher electric fields. 

We now turn from Uni(tOc) to (Tthg(wc). We find 
that saturation can significantly affect THG, as shown 
in Figs. 7(a) and 7(b). Here again the regimes (i) and 
(ii) identified above are relevant. For the results shown 
in Fig. 7(a) we are in regime (i), where both one- and 
two-photon absorption are present. Here both the real 
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FIG. 7. (Color online) Electric field dependence of the non¬ 
linear conductivities o'THG(‘i^c) for different (/i, A), (a) (0,0) 
(squares), (0.2,0) eV (circles), and (0,0.2) eV (diamonds), 
(b) (0.4,0) (squares) and (0,0.4) eV (diamonds). The real 
and imaginary parts are given by filled and hollow symbols, 
and their scales are at the left and right y axis, respectively. 
The dashed curves are drawn to guide the eye. 


and imaginary parts of crTHG(wc) depend strongly on 
the electric field. The imaginary part even changes 
its sign from positive to negative values with increas¬ 
ing the electric field, while the real part shows peaks 
around Eq = 5 x 10^ V/m. Compared to the values 
at Eq = 10® V/m, these peak absolute values are about 
5 times larger. In the measurement of THG®®“^^ the 
procedure used to prepare the samples would indicate 
the chemical potentials should be very low; thus satura¬ 
tion may well occur and the effective THG coefficients 
trTHG(<^c) may be above their perturbative values. For 
the results shown in Fig. 7(b) we are in regime (ii), where 
one-photon absorption is absent but two-photon absorp¬ 
tion is still present. Here the real parts of the crxHG(wc) 
weakly depend on the electric field; however, their imag¬ 
inary parts still strongly depend on the electric field and 
change the sign at about Eq = 10® V/m. 

To qualitatively understand how the saturation affects 
THG, we construct a function 

5(e) = cue; |e|), (24) 

Here Wc, Wc; |/r|) is the analytic perturbative 

third order conductivity of DG®^ at zero temperature, 
with the chemical potential dependence explicitly shown; 
S{e) describes the contribution of the electron states at 
energy e to the THG. For our calculation parameters, 
fiwc = 0.6 eV, T = 300 K, and F^ = Fg = 33 meV, the 
e dependence of 5(e) is shown in Fig. 8. For a given e, 
J^-s S{E)dE is the contribution to the THG of the elec¬ 
trons distributed in the energy range [e — d, e -I- d]. When 
the saturation is induced by the one-photon absorption, 
the electrons are injected into states with energy around 
tujJcl2 from states with energy around —tujJcl2. The 
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FIG. 8. (Color online) 5(e) at fejc = 0.6 eV, T = 300 K, and 
Fi = Fe = 33 meV. Those two vertical dashed lines are at 
e = 0.3 and 0.6 eV respectively. 


contribution of the population changes to the THG is 
approximately oc [5(Swc/2) — 5(—Swc/2)] i?o i with 
originating from the one-photon injection carrier density. 
Similarly, the carriers injected by two-photon absorption 
contribute oc \S{1tjjJc) — S{—fujJc)]EQ, with E^ originat¬ 
ing from the two-photon injection carrier density. Fig¬ 
ure 8 shows the real parts of these two terms are posi¬ 
tive and negative respectively. Thus they give compet¬ 
ing contributions. For the results in Fig. 7(a), at small 
Eq, one-photon absorption dominates, and Re[(TTHG(wc)] 
increases with i?o; at high i?o, two-photon absorption 
starts to play a role, and the appearance of a peak of 
Re[(TTHG(<j-’c)] is possible. The imaginary part and the 
results shown in Fig. 7(b) can also be understood in the 
same way. 


D. Second harmonic generation 

Finally, we consider the dependence of SHG in GG on 
the band gap and electric field amplitude. In parallel 
with our strategy for the third order response, we in¬ 
troduce an effective second-order nonlinear conductivity 
o'shg(wc) which is given by Wc) in the weak 

field limit, and extracted for larger fields from the numer¬ 
ical calculations as sketched in section III. Our results are 
shown in Fig. 9(a) and 9(b). As expected, a nonzero A, 
associated with the lack of centre-of-inversion symmetry, 
leads to a nonzero SHG response. As A is increased from 
0 to 0.4 eV, the real part of cr^^crsHG(wc) decreases from 0 
to a negative minimum value (about —70 pm/V for Eq = 
10® V/m and —30 pm/V for Eq = 2 x 10^ V/m) around 
A = 0.25 eV, then changes sign around A = 0.3 eV 
and reaches a value about 140 pm/V at A = 0.4 eV; 
they show a strong electric field dependence around the 
minimum values. The imaginary part of cr/'^crsHG(wc) 
has positive values with a peak ^ 200 pm/V around 
A = 0.3 eV for both electric field amplitudes consid¬ 
ered. Physically, crsHG(wc) vanishes as A = 0, where the 
centre-of-inversion symmetry is present, and as A —> oo; 
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FIG. 9. (Golor online) (a) A dependence of (tshg(<^c) in GG 
at different electric fields Eq = 10® V/m (squares) and 2 x 
10^ V/m (circles), (b) electric field dependence of crsHG(wc) 
at different A = 0.2 eV (squares) and A = 0.4 eV (circles). 
The real and imaginary parts are given by filled and hollow 
symbols respectively. The dashed lines are drawn to guide 
the eye. 


for large A it vanishes as^^ oc A”'^. Therefore the exis¬ 
tence of a maximum of the magnitude of (Tshg(<Vc) as A 
is increased is not surprising. 

To focus on the electric held dependence of crsHG(wc), 
we plot that dependence in Fig. 9(b) for two gap pa¬ 
rameters, A = 0.2 eV and A = 0.4 eV. The real part 
of crsHG(‘Vc) at A = 0.2 eV shows a strong dependence. 
It changes its sign from negative to positive as the elec¬ 
tric held increases from 0 to 8 x 10^ V/m. The reason 
is similar to the electric held dependence of the third 
order conductivities, and is induced by the saturation ef¬ 
fects. However, the imaginary part changes little over the 
same range of the electric held. For A = 0.4 eV, where 
the saturation effects can be ignored, both the real and 
imaginary parts show minor changes up to a electric held 
20 X 10^ V/m. 

Similar to the estimation for the effective third order 
susceptibilities,®’®^ we calculate the magnitude of the sus¬ 
ceptibility of SHG in GG, by employing the conversion 
of ~ o-shg(wc)/(— 2 /Wceodgr) with the effective 

thickness of graphene dgr = 3.3 A. For maximum values 
of Ict/'^ctshgI 200 pm/V around A = 0.3 eV, we get 
y;G) ^ 2300 pm/V. This value is about 30 times higher 
than the widely used AgGaSe 2 crystal value 68 pm/V at 
the same photon energy,or a few times larger than that 
of monolayer BN, which has a much larger band gap."^®’^"^ 


V. CONCLUSION AND DISCUSSION 

In this work, we numerically solved the semiconduc¬ 
tor Bloch equations, including phenomenological relax¬ 
ation times, for the excitation of both doped and gapped 
graphene excited by a pump pulse, and extracted the 
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effective optical nonlinear conductivities for second har¬ 
monic generation, the Kerr effects, and third harmonic 
generation for a given fundamental photon energy HcOc = 
0.6 eV. We focused on the dependence of these nonlin¬ 
ear coefficients on the chemical potential /i for doped 
graphene, the gap parameter A for gapped graphene, and 
the electric field amplitude for both. We obtained the fol¬ 
lowing results. 

(1) For doped graphene: At weak electric fields, all 
extracted conductivities (both linear and nonlinear) are 
in good agreement with the perturbation results, which 
is a strong evidence of the correctness of both the nu¬ 
merical and perturbation calculations. The numerical 
results also confirm that both the linear dispersion ap¬ 
proximation and the consideration of only optical transi¬ 
tions around the Dirac points are physically appropriate 
in the perturbation calculation with using the standard 
r ■ E interaction^®. With an increase in the electric field 
amplitude, the effective Kerr coefficient shows a depen¬ 
dence on the Held strength, which can be attributed to 
saturation effects. For Hujc > 2|^| where one-photon ab¬ 
sorption exists, the saturation effects can be character¬ 
ized by a saturation field, which for our relaxation param¬ 
eters takes a value of about 3 x 10^ V/m. The amplitude 
of the effective third harmonic generation coefficient can 
increase up to 5 times as the electric field changes from 
10® V/m to 8 X 10^ V/m. However, compared to the 
two orders of magnitude difference between the values 
from the perturbation calculation and experiments,®’®^ 


this small increment indicates that other effects, such as 
the consequences of including more realistic scattering 
and many-body phenomena, may be important. 

(2) For gapped graphene: The third-order optical con¬ 
ductivity for both Kerr effect and third harmonic gener¬ 
ation in gapped graphene shows obvious peaks or valleys 
in its A dependence, which is different from the |/r| de¬ 
pendence in doped graphene due to the nature of the 
velocity matrix elements. The susceptibility of second 
harmonic generation in gapped graphene is of the order 
of 10® pm/V, and shows a complicated dependence on the 
gap parameter A. Compared to the current induced sec¬ 
ond harmonic generation in doped graphene, which could 
be as high as 10"^ pm/V at similar photon energies under 
appropriate conditions,®® the second harmonic genera¬ 
tion coefficients obtained here are smaller but not that 
much. Therefore gapped graphene may also be useful in 
providing a second harmonic generation functionality in 
optical devices. 
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interaction. In our numerical calculation of linear response 
with p ■ A interaction, we find that the inclusion of all k 
in the whole Brillouin zone is necessary for the imaginary 
part of the linear conductivity, even though the inclusion 
of k only with transition energy close to the photon energy 
is adequate for its real part. 



